Delphi Magazine Collection 2001/Delphi Magazine Collection 20001 (2001).iso/DISKS/Issue59/ALFRESCO/AAWidWrd.pas

Pascal/Delphi Source File | 2000-06-05 | 19.6 KB | 747 lines

{*********************************************************}
{* AAWidWrd                                              *}
{* Copyright (c) Julian M Bucknall 2000                  *}
{* All rights reserved.                                  *}
{*********************************************************}
{* Algorithms Alfresco: Wide Word arithmetic             *}
{*********************************************************}

{Note: this unit is released as freeware. In other words, you are free
       to use this unit in your own applications, however I retain all
       copyright to the code. JMB}

unit AAWidWrd;

interface

uses
  SysUtils;

type
  TaaWideWord = class
    private
      FDigitCount : integer;
      FDigits     : PByteArray;
      FPrimeTestCount : integer;
    protected
      procedure wwDivideByDigit(aDigit : byte);
      function wwIsWitness(a, NMinus1 : TaaWideWord) : boolean;
      procedure wwMultiplyByDigit(aDigit : byte);
      function wwNormalize : integer;
      procedure wwRecalcDigitCount;
    public
      constructor Create;
      destructor Destroy; override;

      procedure SetToMax;
      procedure SetToZero;
      function IsDivBy3 : boolean;
      function IsOne : boolean;
      function IsZero : boolean;

      procedure Assign(x : TaaWideWord); overload;
      procedure Assign(x : integer); overload;
      function AsInteger : integer;

      function Compare(x : TaaWideWord) : integer;

      procedure Add(x : TaaWideWord);
      procedure AddOne;
      procedure SubOne;
      procedure Subtract(x : TaaWideWord);
      procedure Multiply(x : TaaWideWord);
      procedure Divide(x : TaaWideWord;
                       aRem : TaaWideWord);

      function IsPrime : boolean;
      function BitCount : integer;
      function GetBit(aInx : integer) : boolean;
      procedure RandomPrime(aBitCount : integer);
      procedure RandomRSAPrime(aBitCount : integer);

      procedure Print;

      property Digits : PByteArray read FDigits;
      property DigitCount : integer read FDigitCount;
      property PrimeTestCount : integer read FPrimeTestCount;
  end;


implementation

{====================================================================}
function MaxI(x, y : integer) : integer;
begin
  if (x < y) then
    Result := y
  else
    Result := x;
end;
{--------}
procedure AddPrim(x, y : PByteArray; aCount : integer);
var
  Temp  : LongWord;
  Carry : LongWord;
  i     : integer;
begin
  Carry := 0;
  for i := 0 to pred(aCount) do begin
    Temp := x^[i] + y^[i] + Carry;
    x^[i] := Temp mod 256;
    Carry := Temp div 256;
  end;
  x^[aCount] := Carry;
end;
{--------}
function ComparePrim(x, y : PByteArray; aCount : integer) : integer;
var
  i : integer;
begin
  for i := pred(aCount) downto 0 do begin
    Result := longint(x^[i]) - y^[i];
    if (Result <> 0) then
      Exit;
  end;
  Result := 0;
end;
{--------}
procedure MultiplyPrim(x, y    : PByteArray;
                   var aXCount : integer;
                       aYCount : integer);
var
  InxX    : integer;
  InxY    : integer;
  MaxX    : integer;
  Carry   : LongWord;
  Temp    : LongWord;
  Result  : array [0..128] of byte;
begin
  {initialize the result}
  FillChar(Result, sizeof(Result), 0);
  {if either x or y had no digits, it was zero; the answer is zero}
  if (aXCount = 0) or (aYCount = 0) then begin
    Move(Result, x^, sizeof(Result));
    aXCount := 0;
    Exit;
  end;
  {for every Y digit we shall multiply by all the X digits}
  MaxX := aXCount;
  for InxY := 0 to pred(aYCount) do begin
    {if the Y digit is zero there'll be nothing to do in this
     loop, so only do work if it is non-zero}
    if (y[InxY] <> 0) then begin
      {start off with a carry of zero}
      Carry := 0;
      {for each X digit, we multiply it with the current Y digit, add
       in the carry from the previous step, and add the current value
       of the result digit}
      for InxX := 0 to pred(MaxX) do begin
        Temp := (LongWord(x[InxX]) * y[InxY]) +
                Result[InxX + InxY] + Carry;
        {strip off the carry and set the result digit}
        Result[InxX + InxY] := Temp mod 256;
        Carry := Temp div 256;
      end;
      {make sure that any remaining carry is saved}
      Result[InxY + MaxX] := Carry;
    end;
  end;
  {return the answer}
  Move(Result, x^, sizeof(Result));
  {return the new number of digits}
  inc(MaxX, aYCount);
  while (x^[pred(MaxX)] = 0) do
    dec(MaxX);
  aXCount := MaxX;
end;
{--------}
procedure PrintPrim(x : PByteArray; aLen : integer);
var
  i, j, k : integer;
begin
  j := 0;
  k := 0;
  for i := 0 to pred(aLen) do begin
    write(Format('%.2x', [x^[i]]));
    inc(j);
    if (j = 4) then begin
      inc(k);
      if (k = 8) then begin
        writeln;
        k := 0;
      end
      else
        write(' ');
      j := 0;
    end;
  end;
  writeln;
end;
{--------}
function SubtractPrim(x, y : PByteArray; aCount : integer) : integer;
var
  Temp   : integer;
  Borrow : integer;
  i      : integer;
begin
  Borrow := 0;
  for i := 0 to pred(aCount) do begin
    Temp := x^[i] - (y^[i] + Borrow);
    if (Temp < 0) then begin
      inc(Temp, 256);
      Borrow := 1;
    end
    else
      Borrow := 0;
    x^[i] := Temp;
  end;
  Result := Borrow;
end;
{====================================================================}


{===TaaWideWord======================================================}
constructor TaaWideWord.Create;
begin
  inherited Create;
  FDigits := AllocMem(129); {enough for 1024 bits plus overflow}
  FPrimeTestCount := 30;    {probability of not being prime: 2^(-30)}
end;
{--------}
destructor TaaWideWord.Destroy;
begin
  if (FDigits <> nil) then
    FreeMem(FDigits);
  inherited Destroy;
end;
{--------}
procedure TaaWideWord.Add(x : TaaWideWord);
begin
  FDigitCount := MaxI(DigitCount, x.DigitCount);
  AddPrim(Digits, x.Digits, DigitCount);
  if (Digits^[DigitCount] <> 0) then
    inc(FDigitCount);
  if (DigitCount > 128) then begin
    SetToMax;
    raise Exception.Create('TaaWideWord.Add: overflow')
  end;
end;
{--------}
procedure TaaWideWord.AddOne;
var
  Temp  : LongWord;
  Carry : LongWord;
  i     : integer;
begin
  Carry := 1;
  i := 0;
  while (Carry = 1) do begin
    Temp := Digits^[i] + Carry;
    Digits^[i] := Temp mod 256;
    Carry := Temp div 256;
    inc(i);
  end;
end;
{--------}
function TaaWideWord.AsInteger : integer;
type
  PInteger = ^integer;
begin
  Result := PInteger(Digits)^;
end;
{--------}
procedure TaaWideWord.Assign(x : TaaWideWord);
begin
  FillChar(Digits^, 129, 0);
  Move(x.Digits^, Digits^, x.DigitCount);
  FDigitCount := x.DigitCount;
end;
{--------}
procedure TaaWideWord.Assign(x : integer);
type
  PInteger = ^integer;
begin
  if (x < 0) then
    raise Exception.Create('TaaWideWord.Assign: cannot convert negative number');
  FillChar(Digits^, DigitCount, 0);
  if (x = 0) then
    FDigitCount := 0
  else begin
    PInteger(Digits)^ := x;
    if (Digits^[3] <> 0) then
      FDigitCount := 4
    else if (Digits^[2] <> 0) then
      FDigitCount := 3
    else if (Digits^[1] <> 0) then
      FDigitCount := 2
    else
      FDigitCount := 1;
  end;
end;
{--------}
function TaaWideWord.BitCount : integer;
var
  Test : byte;
begin
  if (DigitCount = 0) then
    Result := 0
  else begin
    Result := pred(DigitCount) * 8;
    Test := Digits^[pred(DigitCount)];
    if (Test >= $80) then
      inc(Result, 8)
    else if (Test >= $40) then
      inc(Result, 7)
    else if (Test >= $20) then
      inc(Result, 6)
    else if (Test >= $10) then
      inc(Result, 5)
    else if (Test >= $08) then
      inc(Result, 4)
    else if (Test >= $04) then
      inc(Result, 3)
    else if (Test >= $02) then
      inc(Result, 2)
    else
      inc(Result, 1);
  end;
end;
{--------}
function TaaWideWord.Compare(x : TaaWideWord) : integer;
begin
  if (DigitCount < x.DigitCount) then
    Result := -1
  else if (DigitCount > x.DigitCount) then
    Result := 1
  else
    Result := ComparePrim(Digits, x.Digits, DigitCount);
end;
{--------}
procedure TaaWideWord.Divide(x    : TaaWideWord;
                             aRem : TaaWideWord);
var
  Divisor  : TaaWideWord;
  Dividend : TaaWideWord;
  Test     : TaaWideWord;
  TestCompare : integer;
  InxQ     : integer;
  InxX     : integer;
  Factor   : integer;
  SigDigit : byte;
  q        : longint;
begin
  {if the divisor (x) has no digits, it is zero: an error}
  if (x.DigitCount = 0) then
    raise Exception.Create('TaaWideWord.Divide: division by zero');
  {if the divisor (x) is 1, the quotient (us) equals the dividend (us)
   and the remainder is 0}
  if x.IsOne then begin
    aRem.SetToZero;
    Exit;
  end;
  {if the dividend (us) has no digits, it is zero; the quotient (us)
   and remainder must therefore both be zero}
  if (DigitCount = 0) then begin
    aRem.SetToZero;
    Exit;
  end;
  {compare us to x}
  TestCompare := Compare(x);
  {if the dividend (us) is smaller than the divisor (x), the quotient
   (us) is zero and the remainder equals the dividend (us)}
  if (TestCompare < 0) then begin
    aRem.Assign(Self);
    SetToZero;
    Exit;
  end;
  {last special test: if the dividend (us) and divisor (x) are the
   same, the quotient (us) is 1 and the remainder 0}
  if (TestCompare = 0) then begin
    Assign(1);
    aRem.SetToZero;
    Exit;
  end;
  {if we reach this point we actually have to do some work: the
   dividend is greater than the divisor, neither is zero, and the
   divisor is not 1}

  {allocate and initialize some temporaries}
  Test := nil;
  Divisor := nil;
  Dividend := nil;
  try
    Divisor := TaaWideWord.Create;
    Divisor.Assign(x);
    Dividend := TaaWideWord.Create;
    Dividend.Assign(Self);
    Test := TaaWideWord.Create;

    {after the division we shall be holding the quotient; make sure
     we're zeroed out for now}
    SetToZero;

    {make sure the most significant digit of the divisor is greater
     than 128; retain the multiplicative factor to do that}
    Factor := Divisor.wwNormalize;

    {multiply the dividend by the same factor}
    if (Factor <> 1) then
      Dividend.wwMultiplyByDigit(Factor);

    {if the most sigdigit of the dividend is greater than or equal to
     that of the divisor, increment the number of digits in the
     dividend; this'll help once we jump into the division itself}
    if (Dividend.Digits^[pred(Dividend.DigitCount)] >=
        Divisor.Digits^[pred(Divisor.DigitCount)]) then
      inc(Dividend.FDigitCount);

    {Notes:
     let's get some numbers nailed down:
       n = number of digits in divisor
       n+m+1 = number of digits in dividend (the first may be zero)
       m+1 = number of digits in quotient (the first may be zero)
       n = number of digits in remainder
     we calculate the digits in the quotient from most to least
       significant
     the index of the most significant digit in the dividend is given
       by InxX; it starts out as pred(DigitCount) and will reduce by
       one each time through the loop
     the index of the digits in the quotient is given by InxQ, InxQ
       will vary from m down to 0 (ie, m+1 values)
     note that InxQ will be the position of the start of the part of
       the dividend we're looking at; in other words digits InxQ..InxX
       of the dividend form the part of the dividend we're dividing by
       the divisor
    }

    {calculate the position of the most significant digit of both the
     quotient and dividend}
    InxQ := Dividend.DigitCount - Divisor.DigitCount;
    FDigitCount := InxQ;
    InxX := Dividend.DigitCount;
    {get the most significant digit of the divisor}
    SigDigit := Divisor.Digits^[pred(Divisor.DigitCount)];
    {while we are still calculating quotient digits...}
    while InxQ >= 0 do begin
      {calculate our first estimate for this quotient digit q (it may
       be too large of course but the final value will be within 2 of
       this)}
      q := ((LongWord(Dividend.Digits^[InxX]) * 256) +
             Dividend.Digits^[InxX-1]) div SigDigit;
      {refine q if necessary}
      if (q <> 0) then begin
        {it's only a digit so force it in range}
        if (q >= 256) then
          q := 255;
        Test.Assign(Divisor);
        Test.wwMultiplyByDigit(q);
        while (ComparePrim(@Dividend.Digits^[InxQ],
                           Test.Digits,
                           Test.DigitCount+1) < 0) do begin
          dec(q);
          Test.Assign(Divisor);
          Test.wwMultiplyByDigit(q);
        end;
      end;
      {save this digit for the quotient}
      Digits^[InxQ] := q;
      {subtract}
      if (q <> 0) then begin
        SubtractPrim(@Dividend.Digits^[InxQ],
                     Test.Digits,
                     Test.DigitCount+1);
      end;
      {we've done this digit, now do the next}
      dec(InxX);
      dec(InxQ);
    end;
    {we now have the quotient, calculate the number of digits}
    wwRecalcDigitCount;
    {set up the remainder}
    Dividend.wwRecalcDigitCount;
    aRem.Assign(Dividend);
    if (Factor <> 1) then
      aRem.wwDivideByDigit(Factor);
  finally
    Divisor.Free;
    Dividend.Free;
    Test.Free;
  end;
end;
{--------}
function TaaWideWord.GetBit(aInx : integer) : boolean;
const
  Mask : array [0..7] of byte =
           ($01, $02, $04, $08, $10, $20, $40, $80);
var
  ByteNum : integer;
  BitNum : integer;
begin
  ByteNum := aInx div 8;
  BitNum := aInx mod 8;
  Result := (Digits^[ByteNum] and Mask[BitNum]) <> 0;
end;
{--------}
function TaaWideWord.IsDivBy3 : boolean;
var
  i     : integer;
  Count : integer;
  AddIt : boolean;
begin
  Count := 0;
  AddIt := true;
  for i := pred(BitCount) downto 0 do begin
    if GetBit(i) then
      if AddIt then
        inc(Count)
      else
        dec(Count);
    AddIt := not AddIt;
  end;
  Result := (Count mod 3) = 0;
end;
{--------}
function TaaWideWord.IsOne : boolean;
begin
  Result := (DigitCount = 1) and (Digits^[0] = 1);
end;
{--------}
function TaaWideWord.IsPrime : boolean;
var
  a       : TaaWideWord;
  NMinus1 : TaaWideWord;
  TestNum : integer;
  i       : integer;
begin
  {allocate temporaries}
  a := nil;
  NMinus1 := nil;
  try
    a := TaaWideWord.Create;
    NMinus1 := TaaWideWord.Create;
    NMinus1.Assign(Self);
    NMinus1.SubOne;
    {test PrimeTestCount times}
    for TestNum := 1 to PrimeTestCount do begin
      {set a to a random number between 1 and ourselves}
      repeat
        a.SetToZero;
        a.FDigitCount := Random(DigitCount) + 1;
        for i := 0 to pred(a.DigitCount) do
          a.Digits^[i] := Random(256);
        a.wwRecalcDigitCount;
      until (not a.IsZero) and (Compare(a) > 0);
      if wwIsWitness(a, NMinus1) then begin
        Result := false;
        Exit;
      end;
    end;
    Result := true;
  finally
    a.Free;
    NMinus1.Free;
  end;
end;
{--------}
function TaaWideWord.IsZero : boolean;
begin
  Result := DigitCount = 0;
end;
{--------}
procedure TaaWideWord.Multiply(x : TaaWideWord);
begin
  if ((DigitCount + x.DigitCount) > 128) then begin
    SetToMax;
    raise Exception.Create('TaaWideWord.Multiply: overflow');
  end;
  MultiplyPrim(Digits, x.Digits, FDigitCount, x.DigitCount);
end;
{--------}
procedure TaaWideWord.Print;
begin
  if IsZero then
    writeln('zero')
  else
    PrintPrim(Digits, DigitCount);
end;
{--------}
procedure TaaWideWord.RandomPrime(aBitCount : integer);
var
  ByteCount : integer;
  i         : integer;
begin
  {fill with random digits}
  ByteCount := (aBitCount + 7) div 8;
  for i := pred(ByteCount) downto 0 do
    Digits^[i] := Random(256);
  FDigitCount := ByteCount;
  {set the most and least significant bits}
  Digits^[0] := Digits^[0] or $01;
  Digits^[pred(ByteCount)] := Digits^[pred(ByteCount)] or $80;
  {test for being a prime, if not continue adding 2 until we get one}
  while not IsPrime do begin
    AddOne;
    AddOne;
  end;
end;
{--------}
procedure TaaWideWord.RandomRSAPrime(aBitCount : integer);
begin
  repeat
    RandomPrime(aBitCount);
    SubOne;
  until not IsDivBy3;
  AddOne;
end;
{--------}
procedure TaaWideWord.SetToMax;
begin
  FillChar(Digits^, 128, $FF);
  Digits^[128] := 0;
  FDigitCount := 128;
end;
{--------}
procedure TaaWideWord.SetToZero;
begin
  FillChar(Digits^, 129, 0);
  Digits^[128] := 0;
  FDigitCount := 0;
end;
{--------}
procedure TaaWideWord.SubOne;
var
  Done : boolean;
  i    : integer;
begin
  Done := false;
  i := 0;
  while not Done do begin
    if (Digits^[i] = 0) then begin
      Digits^[i] := 255;
      inc(i);
    end
    else begin
      dec(Digits^[i]);
      Done := true;
    end;
  end;
end;
{--------}
procedure TaaWideWord.Subtract(x : TaaWideWord);
var
  Borrow : integer;
begin
  if (DigitCount < x.DigitCount) then
    Borrow := 1
  else
    Borrow := SubtractPrim(Digits, x.Digits, DigitCount);
  if (Borrow <> 0) then
    raise Exception.Create('TaaWideWord.Subtract: negative result');
  wwRecalcDigitCount;
end;
{--------}
procedure TaaWideWord.wwDivideByDigit(aDigit : byte);
var
  Remainder : LongWord;
  Temp      : LongWord;
  i         : integer;
begin
  Remainder := 0;
  for i := pred(DigitCount) downto 0 do begin
    Temp := (Remainder * 256) + Digits^[i];
    Digits^[i] := Temp div aDigit;
    Remainder := Temp mod aDigit;
  end;
end;
{--------}
function TaaWideWord.wwIsWitness(a, NMinus1 : TaaWideWord) : boolean;
var
  d : TaaWideWord;
  Rem : TaaWideWord;
  x : TaaWideWord;
  i : integer;
begin
  {allocate temporaries}
  d := nil;
  x := nil;
  Rem := nil;
  try
    d := TaaWideWord.Create;
    x := TaaWideWord.Create;
    Rem := TaaWideWord.Create;
    d.Assign(1);
    {we're going to be calculating a^(n-1) mod n by the square and
     multiply method}
    for i := pred(NMinus1.BitCount) downto 0 do begin
      x.Assign(d);
      d.Multiply(d);
      d.Divide(Self, Rem);
      if (Rem.DigitCount > Self.DigitCount) then
        writeln('error');
      d.Assign(Rem);
      if d.IsOne then begin
        if (not x.IsOne) and (x.Compare(NMinus1) <> 0) then begin
          Result := true;
          Exit;
        end;
      end;
      if NMinus1.GetBit(i) then begin
        d.Multiply(a);
        d.Divide(Self, Rem);
        if (Rem.DigitCount > Self.DigitCount) then
          writeln('error');
        d.Assign(Rem);
      end;
    end;
    Result := not d.IsOne;
  finally
    d.Free;
    x.Free;
    Rem.Free;
  end;
end;
{--------}
procedure TaaWideWord.wwMultiplyByDigit(aDigit : byte);
var
  i    : integer;
  Temp : LongWord;
  Carry: LongWord;
begin
  if (aDigit <> 1) then begin
    Carry := 0;
    for i := 0 to pred(DigitCount) do begin
      Temp := (LongWord(Digits^[i]) * aDigit) + Carry;
      Digits^[i] := Temp mod 256;
      Carry := Temp div 256;
    end;
    if (Carry <> 0) then begin
      Digits^[DigitCount] := Carry;
      inc(FDigitCount);
    end;
  end;
end;
{--------}
function TaaWideWord.wwNormalize : integer;
var
  Test     : byte;
begin
  Test := Digits^[pred(DigitCount)];
  Result := 1;
  while (Test < 128) do begin
    Result := Result * 2;
    Test := Test * 2;
  end;
  if Result <> 1 then
    wwMultiplyByDigit(Result);
end;
{--------}
procedure TaaWideWord.wwRecalcDigitCount;
begin
  while (DigitCount > 0) and (Digits^[pred(DigitCount)] = 0) do
    dec(FDigitCount);
end;
{====================================================================}

end.